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Abstract — ^ 

The described method solves the inverse problem for supercritical 
blade-to-blade flow on stream surfaces of revolution with variable 
radius and variable stream surface thickness in a relative system. 
Some aspects of shockless design and of leading edge resolution in 
the numerical procedure are depicted. Some supercritical compres- 
sor cascades were designed and their complete flow field results 
were compared with computations of two different analysis methods. 

Nomenclature 



A 

area 

Subscripts : 

BN 

blade number 



La 

Laval number 

1 

upstream 

R 

radius of stream surface 

2 

downstream 

W 

magnitude of velocity vector 

BW 

blade wake 

d 

profile thickness 

SW 

side wall 

h 

stream surface thickness 

ax 

axial 

1 

profile chord length 

m 

meridional 

m 

meridional coordinate 

u 

circumferential 

s 

arc length 

Urn 

Transition 

t 

cascade pitch 



w 

velocity vector 



X 

chord coordinate 



z 

axial coordinate 



r 

circulation 



0 

axial velocity density ratio 



e 

circumferential angle 



>3 

flow angle 



€ 

inclination angle of stream surface 



P 

density 



<P 

potential function 



U> 

stream function 



C 0 

angular velocity 
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INTRODUCTION 


Increasing requirements on turbomachines concerning efficiency, 
compact construction and density of power lead to aerodynamically 
highly loaded blades. The admissible blade load and the profile 
losses are determined by the boundary layer development. High 
pressure ratios per stage and high turning of the flow increase 
the risk of boundary layer separation with the result of strongly 
growing losses. This problem is intensified by the risk of arising 
compression shocks. They are caused by the supercritical through- 
flow (with local supersonic regions) , which is necessary for high 
mass flow density. 

In regions with pressure rise, boundary layer separation can only 
be avoided by careful blade profiling for flow with minimum loss. 
In the past , turbomachinery bladings have mostly been designed 
with the aid of profile families. But in this way, depending on 
the plurality of parameters, a loss minimization is not possible. 
Especially in the transonic velocity region this procedure is in- 
sufficient since shockfree solutions can be found with only poor 
chances by iterative contour variation. In this region of maximum 
mass flow density very small variations of the geometry are con- 
nected with very high changes in the flow velocity. Therefore, it 
is convenient to prescribe the physical quantity, where the great 
changes appear , and to calculate the small but important varia- 
tions of the geometry. 

This alternative is given by inverse design: Starting from a pre- 
scribed shockfree velocity distribution, the corresponding profile 
contour is calculated numerically. By this means a perfect tai- 
loring of the blade to the required turning problem is possible. 

Up to now a perfect three-dimensional inverse design of flow 
fields in turbomachines is not executed since this problem is 
quite overdetermined. At present a standard procedure is to start 
a quasi three-dimensional computation by calculating the flow on 
meridional planes by an analysis code (duct- or through-flow) 
to get the starting values of the calculation on several blade-to- 
blade planes (S^ distributed along the blade height. 

This multi-section design of the blade can be realized by inverse 
computation. The following inverse computation method is an ex- 
tension of the former cylindrical version [1,2] to the design on 
stream surfaces of revolution with variable radius and variable 
stream surface thickness in a relative system. This development is 
a further step to approximate the real physical behaviour of the 
flow. The method is applied to the multi-section design of a 
three-stage research compressor which is now in construction. 
Computations for comparison were carried out with two different 
analysis codes. 
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FUNDAMENTAL EQUATIONS OF THE METHOD 

Since inverse design strives for low loss flow without shocks and 
boundary layer separation, Prandtl's concept of distinct potential 
flow and boundary layer calculation is applied. In Fig. 1 the 
fundamental process of the method is sketched. The method computes 
the steady compressible potential flow in the passage between two 
blades of unknown shape from far upstream to far downstream on a 
stream surface of revolution. Besides the upstream and downstream 
velocity vector, the velocity distributions are prescribed along 
the arc length of the stagnation streamlines with their periodic 
parts in the upstream and downstream regions and along the blade 
suction and pressure sides. In this way velocity gradients can be 
prescribed, which is important for the boundary layer development 
and a prerequisite for loss minimization. Moreover, radius and 
thickness of the stream surface of revolution are prescribed along 
the axial coordinate. These boundary conditions are transformed by 
integration into the computation plane with stream function 
coordinates and their normals. The computation grid is rectangular 
in this plane and contour adapted in the physical plane. 
Therefore, no interpolations are necessary on the boundaries. 

The equations of continuity and motion for steady, isentropic flow 
on a stream surface of revolution are 



d9 dm dm 


with the isentropic relation 

p = f (W, u, R) (3) 

The appearance of a variable stream surface radius in the 
fundamental equations requires a different treatment of rotor and 
stator flows since in the relative system of the rotor an energy 
alteration is connected with a radius alteration. So the relative 
velocity W can no longer be prescribed by a potential. In the ab- 
solute system {following Vavra [3]) an equivalent potential <p can 
be defined by A 

7 • <P A = = w + i^a> R and 7 * v^ = 0 . ( 4 ) 

But the contour velocity distribution has to be prescribed in the 
relative system, thus one coordinate direction is given by the 
streamlines in the relative system. Together with the potential 
lines normal to the absolute velocity an oblique-angled coordinate 
system results, Fig. 2. 
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At velocities W ■■ wR and flow angles ft * n/2 (e.g. in the stagna- 
tion point region) the angle a approaches zero so that both co- 
ordinate directions coincide. Following from numerical reasons 
this system consisting of an absolute potential and relative 
stream function is inconvenient for use as a computational grid. 

The potential-streamf unction-plane is the computation plane of the 
inverse design method. Furthermore, for using it in the rotor the 
equation of motion is reduced to 7 > w = 0. The consideration of 
the rotational character of the flow occurs by variation of the 
total quantities dependent on the stream surface radius. (Another 
consideration by definition of a transformed potential is pub- 
lished in [4,5].) 

Since the critical sonic velocity is no longer constant, because 
of the variation of the total temperature, it is no longer appli- 
cable for normalization of the velocity like in the stator case. 
Hence, the upstream velocity W. is now applied for this purpose. 
The decision which difference operator for consideration of the 
type-dependence of the differential equation system has to be 
applied is taken by the magnitude of the local Mach number. 

By transformation of the fundamental equations into the potential- 
streamf unction-plane and elimination of the flow angle the full 
potential equation follows [6] : 


C 


1 


o * 

d In W 


' <P 


+ C. 


2 * 
® In w 

din W*1 

2 

4 . C 

f din W*1 

.2 + C 3 

o v 

8 <P J 

C 4 

l 9 v J 


with C 1 . . . C 7 
and W 


+ C, 


din W 

d <p 


+ C, 


din W 

dt 


+ C. 


0 


f (W , u), R , h) (5) 

W/W^ for the rotor 

La for the stator 


The flow field is computed by the solution of the corresponding 
difference equation system applying relaxation combined with 
multi-grid. The change of type (elliptic-hyperbolic) from sub- 
sonic to supersonic flow regions depending on the sign of the 
coefficient is considered by modified difference equations. The 
transformation of the solution back into the physical plane is 
performed by integration of the equations of continuity and 
motion. It yields the field boundaries, i.e. the blade profiles 
and the cascade geometry. 


For solving eq.(5) the dependence of the radius on the potential 
and stream function R = f{<p,w ) is necessary. Since only the axial 
development R = f(z) is known by the prescription, this relation 
can only be discovered by an additional iteration in the course of 
the solution process. The radius distribution has to fulfill the 
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condition of constant values in circumferential direction. Normal- 
ly the same is true for the stream surface thickness. Moreover, 
the desired values of turning angle, pitch-chord ratio or blade 
thickness distributions are attainable by iterative variations of 
the prescribed boundary values (which are selfacting included in 
the code) . The whole geometry of the problem is always the result 
of the computation and therefore completely unknown at the 
beginning . 


CIRCULATION, LEADING EDGE RESOLUTION, SHOCKLESS DESIGN 


The profile circulation necessary for the desired turning, follows 
from ( s . Fig . 3 ) : 


r 


p 


w • ds 


■ cos /? ^ 


2 nrR^ 

BN 


W- ■ cos 13 


2nR 2 

BN ' Q BW 


( 6 ) 


The line integral of the velocity along the computation grid 
boundary can be converted by the law of Stokes into an area inte- 
gral 

^ w • ds = - 2 | u 1 dA . ( 7 ) 

The rotational vector w indicates in axial direction, dA is 
perpendicular to the through-flowed area_. For^ the stator (w = 0) 
and for a rotor with constant radius (« 1 dA) the value of the 
line integral equals zero. Because of the reduction of the equa- 
tion of motion this is also true for rotor flow with varying 
radius. The additionally existing circulation inside the computa- 
tion grid is thereby neglected. 


w • ds 


0 


( 8 ) 


The profile circulation which is necessary for the actual turning 
problem and which should be rendered by the prescribed velocity 
distribution at the beginning of the design process is 


r 


p 


2ti- R^ 

“bn 




W, 


• cos/? 




• cos /?^ 


(9) 


The circulation inside the computation grid neglected in the rotor 
case with varying radius can be estimated in maximum if the flow 
conus area divided by the blade number (i.e. vanishing profile 
area) is assumed as upper limit for the integration area: 


2 

n n 
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w • ds 


BN 


( 10 ) 
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The relative deviation (refered to the circulation 
equation of motion) is 


AT 


- 2 - 


W. 


-cos B. 


n ■ n • 


R 
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COS Vw 


2 

1 



of 


the reduced 


( 11 ) 


In the case of compressor cascades ( B 1 2 90 ) the reduced equation 
of motion yields lower circulation for increasing stream surface 
(R 2 /Ri > 1) and higher circulation for decreasing stream surface. 
In standard cases the deviation amounts to less than 5 percent 
according to a turning angle deviation of less than 1 degree. 


The leading edge region of a profile has special requirements for 
the numerical aspects of a computer program for calculation of the 
flow around an airfoil. In the design method this difficulty be- 
comes especially clear since even the prescription data - the 
velocity distribution on the boundaries of the flow field to be 
computed - show the strong gradients in the stagnation point 
region (Fig. 4). This area can be recorded only insufficiently in 
an equidistant divided computation grid. 


For appropriate resolution of the blade nose region it was found 
that the number of points on the flow field boundaries should be 
up to 2 or 2 times higher than that of the normal grid. Thereby 
local grid refinement is provided for the regions with steep gra- 
dients. For smaller point distances the possibility of emboxing of 
refinements was established. In a corresponding fitted arrangement 
a gradual transition of the mesh size follows. This is especially 
favourable for the accuracy of the solution. For even higher accu- 
racy, a feedback calculation can be performed which uses the 
results of the fine grid for recalculation in the coarse grid in 
an iterative way with overlapping boundaries of both regions. 

In case of velocity prescriptions on the boundaries of local 
supersonic regions an "ill-posed problem" is treated, i.e. no 
physical solution may exist. Numerically this often leads to the 
formation of oscillating shocks in the flow field, shown in 
Fig. 5. If they are weak enough, a provision for cancellation of 
these shocks is given in Fig. 6: Following the plotted character- 
istic directions in the supersonic region from the concerned 
region to the corresponding boundary values, these values can be 
modified for generation of additional expansion waves to remove 
the shocks [7] . 


RESULTS 


The first example is a cascade for a compressor stator hub section 
with supercritical flow (local supersonic region on the suction 
side), Fig. 7. The high subsonic velocity La. = 0.90 is deceler- 
ated to the downstream value La_ = 0.593 by arurning angle of 25 

degrees and a relatively high pitch-chord ratio of t/1 =1.0 . 

ax 
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The stream surface radius R increases by 30 percent from upstream 
to downstream and the stream surface thickness h decreases by 25 
percent in the same direction. Their prescribed slopes dependent 
on the axial length z are given in Fig. 8. The curves of the inner 
and outer radius of the stream surface of revolution consist of 
cosine slopes, the maximum angle of inclination of the stream sur- 
face is c = 25 degrees. In the cascade region the stream surface 
thickness distribution follows a cubic parabola, in the upstream 
and downstream region it is calculated by constant flow area. 

In the upper part of Fig. 7 the full line shows the prescribed 
velocity distribution on the blade suction and pressure side. This 
roof top distribution with maximum Mach number of 1.19 (La = 1.15) 
was chosen for separation-free flow with high loading. On the suc- 
tion side transition takes place at the beginning & of the pressure 
rise at 31 percent of chord length (at Re = 4.7 10 and Tu = 4%) . 

In the lower part of Fig. 7 the computed profile shape is plotted. 
The dashed line marks the contour of the potential flow computa- 
tion from which the manufacturing contour (full line) is derived 
by subtraction of the boundary layer displacement thickness (com- 
puted by Rotta ' s integral method [8]). The complete cascade geome- 
try and the flow field characterized by the (full) lines of con- 
stant velocity (with the sonic line La = 1.0) are shown in Fig. 9. 
The cascade geometry data were used as input for the analysis code 
of Lucking [9]. The results, the contour velocity distribution 
(crosses in Fig. 7) and the velocity distributions in the flow 
field (dashed lines in Fig. 9) agree well with the distributions 
of the inverse code, even in the supersonic region. 

Moreover, in Fig. 10 the course of the lines of constant radius 
(dashed-dotted), which can only be calculated iteratively (see 
above), coincide well with the demanded circumferential direction. 

In Fig. 11 the prescribed velocity distribution of a rotor tip 
section (n = 3600 rpm) is plotted together with the resulting 
profile shape (both full lines). Despite of the low turning of 6 
degrees, due to the high upstream velocity of La = 0.866 and the 
high pitch-chord ratio of t/1 = 1,761 the loading is high enough 
to require local supersonic flow on the suction side. In Fig. 12 
the cascade geometry and the flow field consisting of lines of 
constant velocity is drawn, showing the great stagger of this 
design. The geometry of this result was again used to compute the 
velocity distributions for comparison by the finite volume method 
originating in P.W. McDonald [10] . The crosses in Fig. 11 and the 
dashed lines in Fig. 12 exhibit satisfying agreement with the 
design computation (full lines). 

The velocity distribution of the rotor hub section (Fig. 13) 
belonging to the preceding rotor tip section was mainly influenced 
by a desired maximum thickness of the profile (in consideration of 
structural reasons). Therefore, high velocities appear on the 
suction and pressure side without high aerodynamic loading . The 
design on stream surface with radius increase of 14 percent (full 
line) is compared with plane flow design for equal upstream and 
downstream velocity vectors. It is to be seen that in the plane 
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flow case higher circulation is needed for the same turning prob- 
lem but a thinner profile results compared to the design on 
increasing radius. In Fig. 14 the cascade geometry is demonstrated 
and the isolines of the velocity are compared with the results of 
the finite-volume method [10] . In the front part of the flow 

channel including the local supersonic region the velocity field 
compares well. In the rear part, refering to a local aft-accelera- 
tion behind the supersonic patch the coincidence of the isolines 
is somewhat disturbed. 

CONCLUSION 

The present extended inverse method seems to be an effective 

procedure to design highly loaded axial compressor cascades on 
stream surfaces of revolution. It produces accurate results 
compared with complete flow field results of other methods and was 
successfully applied to cascade and multi-section compressor blade 
design . 
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Fig. 7 : 

Comparison of Laval number distributions on the 
blade between design and analysis calculation on 
stream surfaces of revolution for the stator hub 
section. Additionally the potential flow contour 
{dashed line) and the metal section contour 
(full line) are indicated. 


Slop** of inn*r and outer radius of the stream 
surface (top) and stream surface thickness 
(bottom) for the stator hub section. 




* j distributions of Laval numbers between 

on stream surfaces of revolution for the stator hub 
section (increment ALa * 0.05, sonic line L 
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Fiq. 11 : 

Comparison of Laval number distributions on the 
blade between design and analysis calculation on 
stream surfaces of revolution for the rotor tip 
section . 






calculation (full liaa, aonic line daahed- 
dottad) and analysis calculation (dashed line, 
■onic line abort daahed) for the rotor tip 
aection (increment ALa - 0.51. 
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